function data = PV_dgcomp_series(fname,itime,N,k,l,outfile,varargin)
% data = PV_dgcomp_series(fname,itime,N,k,l,outfile)
% data = PV_dgcomp_series(fname,itime,N,k,l,outfile,select)
%
% Call PV_dgcomp for a series of files and prepare a .pvd file to
% collect the generated vtk files.
%
% fname nust contain a '%' to indicate the subdomain number and a '$'
% to indicate the suffix 'grid', 'base' or 'res-XXXX'.
%
% The optional argument 'select' can be used as follows:
% C|P : conservative or primitive variables
% D|T : deviations or total values
%
% Moreover, the collected and/or postprocessed data are collected in
% the output argument data in a format which can be passed to
% dgcomp_interpolation to interpolate at arbitrary points.

 part_format = '%.4i'

 % input arguments
 if(nargin>6)
   select = varargin{1};
 else
   select = '';
 end
 primitive = not(isempty(strfind(select,'P')));
 total     = not(isempty(strfind(select,'T')));

 % Summary file
 fnameNX = [outfile,'.pvd'];
 fid = fopen(fnameNX,'w');
 
 line = '<?xml version="1.0"?>';
 fprintf(fid,'%s\n',line);
 line = '<VTKFile type="Collection" version="0.1" byte_order="BigEndian">';
 fprintf(fid,'%s\n',line);
 line = ' <Collection>';
 fprintf(fid,'%s\n',line);

 % Part loop
 Npos = find(fname=='%');
 for j=1:length(itime)
   for i=0:N
     fnameN = [fname(1:Npos-1),num2str(i,part_format),fname(Npos+1:end)];
     Xpos = find(fnameN=='$');
     fnameNX = [fnameN(1:Xpos-1),'grid',fnameN(Xpos+1:end)];
     grid = load(fnameNX);
     fnameNX = [fnameN(1:Xpos-1),'base',fnameN(Xpos+1:end)];
     base = load(fnameNX);
     fnameNX = [fnameN(1:Xpos-1),'res-',num2str(itime(j),'%.4i'), ...
                fnameN(Xpos+1:end)];
     res = load(fnameNX);

     if(primitive)
       varnames = {'pressure','temperature','velocity'};
       for it=1:size(res.uuu,1)-2-(grid.grid.d)
         varnames{end+1} = ['tracer ',num2str(it)];
       end
       varnames{end+1} = 'pi';
       varnames{end+1} = 'theta';
       %for it=1:size(res.add_diags,1)
       %  varnames{end+1} = ['diag ',num2str(it)];
       %end
     else
       varnames = {'density','energy','momentum'};
       for it=1:size(res.uuu,1)-2-(grid.grid.d)
         varnames{end+1} = ['tracer ',num2str(it)];
       end
       %for it=1:size(res.add_diags,1)
       %  varnames{end+1} = ['diag ',num2str(it)];
       %end
     end

     if(or(total,primitive)) % need some postprocessing
       fnameNX = [fnameN(1:Xpos-1),'init',fnameN(Xpos+1:end)];
       ref = load(fnameNX);
     end

     if(total==1)
       if(primitive==1)
         [tmp,tmp,uuu,pitheta] = dgcomp_postp(res.uuu , base.base , ...
                                ref.uuu_ioref , ref.atm_ref , ref.phc);
	 %diags = [pitheta ; res.add_diags];
       else
         uuu = res.uuu + ref.uuu_ioref;
         uuu(1,:,:) = uuu(1,:,:) + permute(ref.atm_ref.rho,[3 1 2]);
         uuu(2,:,:) = uuu(2,:,:) + permute(ref.atm_ref.e,[3 1 2]);
	 %diags = res.add_diags;
       end
     else
       if(primitive==1)
         [tmp,uuu,tmp,pitheta] = dgcomp_postp(res.uuu , base.base , ...
                                ref.uuu_ioref , ref.atm_ref , ref.phc);
	 % This is a trick to compute pi and theta for the reference
	 % state, so that we can plot the deviations.
         [tmp,tmp,tmp,io_pith] = dgcomp_postp(0*res.uuu,base.base , ...
                                ref.uuu_ioref , ref.atm_ref , ref.phc);
	 %diags = [pitheta-io_pith ; res.add_diags];
       else
         uuu = res.uuu;
	 %diags = res.add_diags;
       end
     end
                  diags = zeros(0,size(uuu,2),size(uuu,3)); % change this!!!

     % Add the Courant numbers
     diags(end+1,1,:) = 1/base.base.p(1,1)*permute(res.c_tot,[3 1 2]);
     diags(end+1,1,:) = 1/base.base.p(1,1)*permute(res.c_adv(1,:),[3 1 2]);
     diags(end+1,1,:) = 1/base.base.p(1,1)*permute(res.c_adv(2,:),[3 1 2]);
     % Following modes are left to 0
     varnames{end+1} = 'Courant (tot)';
     varnames{end+1} = 'Courant (adv 1)';
     varnames{end+1} = 'Courant (adv 2)';

     fnameNX = [outfile,'-',num2str(itime(j)), ...
                        '-',num2str(i,part_format),'.vtu'];
     PV_dgcomp(grid.grid,base.base,uuu,k,l,fnameNX,res.test_name, ...
               diags,varnames);

     line = ['  <DataSet timestep="', num2str(res.time), ...
             '" group="" part="', num2str(i,'%i'), ...
             '" file="',fnameNX,'"/>'];
     fprintf(fid,'%s\n',line);

     if(nargout>0) % collect data
       data(j).grid{i+1}     = grid.grid;
       data(j).ddc_grid{i+1} = grid.ddc_grid;
       data(j).base{i+1}     = base.base;
       data(j).uuu{i+1}      = [uuu;diags];
     end
   end
 end

 line = ' </Collection>';
 fprintf(fid,'%s\n',line);
 line = '</VTKFile>';
 fprintf(fid,'%s\n',line);

 fclose(fid);

return

